import sys
sys.path = ['/home/as259691/PycharmProjects/FluidDyn1D'] + sys.path
from src.main import *
from src.plot_fields import *
%matplotlib inline
rc('figure', figsize=(10,5))
rc('figure', dpi=100)
Ici on va réaliser une simulation sans diffusion pour différentes écritures de notre équation thermique.
La résolution se fait à chaque fois en WENO avec Euler explicite en temps.
# d = 6./100*Delta/2.
phy_prop = PhysicalProperties(Delta=0.02, v=0.2, dS=0.005**2,
lda1=5.5*10**-2, lda2=15.5, rho_cp1=70278., rho_cp2=702780., diff=1.,
alpha=0.06, a_i=357.)
markers = Bulles(phy_prop=phy_prop)
Formulation = [Problem, ProblemConserv2]
n = 1000
t_fin = 0.2
def compare_energy_forme(formu, phy_prop, num_prop, markers, t_fin):
fig1,ax1 = plt.subplots(1)
ax1.set_title('Énergie en fonction du temps')
for form in formu:
print()
prob = form(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob.energy
print(prob.name)
print('==========================')
t, e = prob.timestep(t_fin=t_fin, number_of_plots=5, debug=False, plotter=Plotter('decale'))
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob.name)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob.dt / E0 # on a mult
# par Dt / rho_cp_l T_l V
print('dE*/dt* = %f' % dedt_adim)
le = fig1.legend()
En fait s'il n'y a pas de convection il n'y a pas de différence entre les différentes formes, à l'exception de la moyenne utilisée pour $\frac{1}{\rho C_p}$
num_prop = NumericalProperties(dx=3.9*10**-5, schema='weno', time_scheme='rk4', phy_prop=phy_prop, cfl=0.5)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 6.918433404737903e-06 Cas : mixte, rk4, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = -0.000006 dt fourier 6.918433404737903e-06 Forme conservative boniou, Cas : mixte, rk4, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = 0.000022
num_prop = NumericalProperties(dx=3.9*10**-5, schema='weno', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 6.918433404737903e-06 Cas : mixte, euler, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = -0.000005 dt fourier 6.918433404737903e-06 Forme conservative boniou, Cas : mixte, euler, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = 0.000023
num_prop = NumericalProperties(dx=3.9*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 6.918433404737903e-06 Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = -0.000006 dt fourier 6.918433404737903e-06 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = 0.000005
dx_list = [3.9*10**-6, 3.9*10**-5, 1*10**-5]
t_fin = 0.01
for dx in dx_list:
num_prop = NumericalProperties(dx=dx, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 6.896863867106643e-08 Cas : mixte, euler, weno upwind, dx = 3.90016e-06, dt = 6.89686e-08, cfl = 0.00353671 ========================== dE*/dt* = -0.000000 dt fourier 6.896863867106643e-08 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 3.90016e-06, dt = 6.89686e-08, cfl = 0.00353671 ========================== dE*/dt* = 0.000000 dt fourier 6.918433404737903e-06 Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = -0.000103 dt fourier 6.918433404737903e-06 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224 ========================== dE*/dt* = 0.000034 dt fourier 4.538601983461999e-07 Cas : mixte, euler, weno upwind, dx = 1.0005e-05, dt = 4.5386e-07, cfl = 0.00907267 ========================== dE*/dt* = -0.000004 dt fourier 4.538601983461999e-07 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 1.0005e-05, dt = 4.5386e-07, cfl = 0.00907267 ========================== dE*/dt* = 0.000001
dx_list = [7*10**-5]
t_fin = 0.01
for dx in dx_list:
num_prop = NumericalProperties(dx=dx, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 2.232841866976439e-05 Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05, cfl = 0.063636 ========================== dE*/dt* = -0.000469 dt fourier 2.232841866976439e-05 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05, cfl = 0.063636 ========================== dE*/dt* = 0.000172
Ici on ne change pas $\Delta x$, mais on diminue dt_min pour qu'il soit contraignant
dt_list = [2*10**-5, 1*10**-5, 5*10**-6]
t_fin = 0.01
for dt in dt_list:
num_prop = NumericalProperties(dx=7*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop, dt=dt)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt min 2e-05 Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2e-05, cfl = 0.057 ========================== dE*/dt* = -0.000419 dt min 2e-05 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2e-05, cfl = 0.057 ========================== dE*/dt* = 0.000154 dt min 1e-05 Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 1e-05, cfl = 0.0285 ========================== dE*/dt* = -0.000209 dt min 1e-05 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 1e-05, cfl = 0.0285 ========================== dE*/dt* = 0.000076 dt min 4.9999999999999996e-06 Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 5e-06, cfl = 0.01425 ========================== dE*/dt* = -0.000104 dt min 4.9999999999999996e-06 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 5e-06, cfl = 0.01425 ========================== dE*/dt* = 0.000038
t_fin = 0.01
phy_prop = PhysicalProperties(Delta=0.02, v=0., dS=0.005**2,
lda1=5.5*10**-2, lda2=15.5, rho_cp1=70278., rho_cp2=702780., diff=1.,
alpha=0.06, a_i=357.)
num_prop = NumericalProperties(dx=7*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 2.232841866976439e-05 Cas : diffusion, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05 ========================== dE*/dt* = 0.000010 dt fourier 2.232841866976439e-05 Forme conservative boniou, Cas : diffusion, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05 ========================== dE*/dt* = -0.000000
t_fin = 0.01
phy_prop = PhysicalProperties(Delta=0.02, v=0.2, dS=0.005**2,
lda1=5.5*10**-2, lda2=15.5, rho_cp1=70278., rho_cp2=702780., diff=0.,
alpha=0.06, a_i=357.)
num_prop = NumericalProperties(dx=7*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)
dt fourier 2.232841866976439e-05 Cas : convection, euler, weno upwind, dx = 7.01754e-05, cfl = 0.063636 ========================== dE*/dt* = -0.000574 dt fourier 2.232841866976439e-05 Forme conservative boniou, Cas : convection, euler, weno upwind, dx = 7.01754e-05, cfl = 0.063636 ========================== dE*/dt* = 0.000166
Schemas = ['upwind', 'center', 'weno', 'weno upwind']
Time_scheme = ['euler', 'rk4']
def compare_energy_schema(schemas, form, time_scheme, phy_prop, markers, t_fin):
fig1,ax1 = plt.subplots(1)
ax1.set_title('Énergie en fonction du temps')
for schem in schemas:
for ts in time_scheme:
a = Plotter('decale')
num_prop = NumericalProperties(dx=3*10**-5, schema=schem, time_scheme=ts, phy_prop=phy_prop)
print()
prob = form(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob.energy
print(prob.name)
print('==========================')
t, e = prob.timestep(t_fin=t_fin, number_of_plots=5, debug=False, plotter=a)
a.ax.set_xlim(0., phy_prop.Delta/2)
a.ax.set_ylim(0.6, 1.1)
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob.name)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob.dt / E0 # on a mult
# par Dt / rho_cp_l T_l V
print('dE*/dt* = %f' % dedt_adim)
le = fig1.legend()
t_fin = 0.01
fig1, ax1 = plt.subplots(1)
a = Plotter('decale')
num_prop = NumericalProperties(dx=1*10**-5, schema='weno', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = Problem(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
a.ax.set_ylim(0.5,1.05)
a.ax.set_xlim(0., phy_prop.Delta/4.)
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob_ref.name)
ax1.legend()
ax1.grid(b=True)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt / E0 # on a mult
num_prop = NumericalProperties(dx=5*10**-6, schema='weno', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = Problem(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
a.ax.set_ylim(0.5,1.05)
a.ax.set_xlim(0., phy_prop.Delta/4.)
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob_ref.name)
ax1.legend()
ax1.grid(b=True)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt / E0 # on a mult
# par Dt / rho_cp_l T_l V
print('dE*/dt* = %f' % dedt_adim)
dt fourier 4.538601983461999e-07 Cas : mixte, rk4, weno, dx = 1.0005e-05, dt = 4.5386e-07, cfl = 0.00907267 ========================== dt fourier 1.1335161290322582e-07 Cas : mixte, rk4, weno, dx = 5e-06, dt = 1.13352e-07, cfl = 0.00453406 ========================== dE*/dt* = -0.000000
t_fin = 0.01
Schemas = ['upwind', 'weno', 'weno upwind']
compare_energy_schema(Schemas, Problem, Time_scheme, phy_prop, markers, t_fin)
dt fourier 4.0888316389624144e-06 Cas : mixte, euler, upwind, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = -0.000083 dt fourier 4.0888316389624144e-06 Cas : mixte, rk4, upwind, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = -0.000084 dt fourier 4.0888316389624144e-06 Cas : mixte, euler, weno, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = -0.000024 dt fourier 4.0888316389624144e-06 Cas : mixte, rk4, weno, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = -0.000026 dt fourier 4.0888316389624144e-06 Cas : mixte, euler, weno upwind, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = -0.000054 dt fourier 4.0888316389624144e-06 Cas : mixte, rk4, weno upwind, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = -0.000055
Schemas = ['upwind', 'weno', 'weno upwind']
Time_scheme = ['euler', 'rk4']
compare_energy_schema(Schemas, ProblemConserv2, Time_scheme, phy_prop, markers, t_fin)
dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : mixte, euler, upwind, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = 0.000033 dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : mixte, rk4, upwind, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = 0.000032 dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : mixte, euler, weno, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = 0.000036 dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : mixte, rk4, weno, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = 0.000035 dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = 0.000018 dt fourier 4.0888316389624144e-06 Forme conservative boniou, Cas : mixte, rk4, weno upwind, dx = 3.003e-05, dt = 4.08883e-06, cfl = 0.0272316 ========================== dE*/dt* = 0.000018
t_fin = 0.2
fig1, ax1 = plt.subplots(1)
for dx in [2*10**-5, 4*10**-5, 7*10**-5, 10*10**-5]:
print('~~~~~~~~~~~~~~~~~~~~~~~~~~')
a = Plotter('decale')
num_prop = NumericalProperties(dx=dx, schema='weno', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = Problem(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy_m
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
l = ax1.plot(t, e, label=prob_ref.name)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt / E0 # on a mult
print('dE/dt = %g' % dedt_adim)
num_prop = NumericalProperties(dx=dx, schema='weno upwind', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = ProblemConserv2(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy_m
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
l = ax1.plot(t, e, label=prob_ref.name)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt
print('dE/dt = %g' % dedt_adim)
a.ax.set_ylim(0.,1.05)
a.ax.set_xlim(0., phy_prop.Delta/4.)
ax1.legend()
ax1.grid(b=True)
~~~~~~~~~~~~~~~~~~~~~~~~~~ dt fourier 7.283608525474247e-06 Cas : mixte, rk4, weno, dx = 4.00802e-05, dt = 7.28361e-06, cfl = 0.0363452 ========================== dE/dt = -5.40692e-12 dt fourier 7.283608525474247e-06 Forme conservative boniou, Cas : mixte, rk4, weno upwind, dx = 4.00802e-05, dt = 7.28361e-06, cfl = 0.0363452 ========================== dE/dt = 1.6338e-07 ~~~~~~~~~~~~~~~~~~~~~~~~~~ dt fourier 2.232841866976439e-05 Cas : mixte, rk4, weno, dx = 7.01754e-05, dt = 2.23284e-05, cfl = 0.063636 ========================== dE/dt = -1.8396e-11 dt fourier 2.232841866976439e-05 Forme conservative boniou, Cas : mixte, rk4, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05, cfl = 0.063636 ========================== dE/dt = 5.1594e-07 ~~~~~~~~~~~~~~~~~~~~~~~~~~ dt fourier 4.534064516129032e-05 Cas : mixte, rk4, weno, dx = 0.0001, dt = 4.53406e-05, cfl = 0.0906813 ========================== dE/dt = -3.69751e-11 dt fourier 4.534064516129032e-05 Forme conservative boniou, Cas : mixte, rk4, weno upwind, dx = 0.0001, dt = 4.53406e-05, cfl = 0.0906813 ========================== dE/dt = 1.16584e-06